using System;
using System.Collections.Generic;
using System.Runtime.Serialization;
using System.Text;


namespace Aima.Core.Util.math
{
/**
 * Matrix class.
 * provides the fundamental operations of numerical linear algebra
 * Various constructors create Matrices from 2D arrays of floats
 * property accessors to submatrices and matrix elements. 
 * methods for basic matrix arithmetic, norms and elemen ops
 * Methods for reading and printing matrices 
 * Complex matrices may be handled in a future version.
 * 
 * Five fundamental matrix decompositions - consist of pairs or triples of
 * matrices, permutation vectors, produce results in five decomposition classes
 * to compute solutions of simultaneous linear equations, determinants, inverses and other matrix functions. The five decompositions are:
 *      Cholesky Decomposition of symmetric, positive definite matrices.
 *      LU Decomposition of rectangular matrices.
 *      QR Decomposition of rectangular matrices.
 *      Singular Value Decomposition of rectangular matrices.
 *      Eigenvalue Decomposition of both symmetric and nonsymmetric square matrices.
 * apparenlty from The MathWorks, Inc. and the National Institute of Standards and Technology. v5 August 1998
 */

    [Serializable]
    public class Matrix : ICloneable, ISerializable
    {
        // internal array storage.
        private readonly double[][] _a;
	    // Row and column dimensions.
        private readonly int _m;
        private readonly int _n;

        public double[][] MatrixArray
        {
            get { return _a; }
        }

        public int RowDimension
        {
           get{return _m;}
        }
     
        public int ColumnDimension
        {
            get{return _n;}
        }


        #region Constructors
        // ctor - Construct an m-by-n matrix of zeros.
        public Matrix(int m, int n)
        {
            _m = m;
            _n = n;
            _a = ReturnRectangularDoubleArray(m, n);
        }


        // ctor - Construct an m-by-n constant matrix - Fill the matrix with this scalar value.
        public Matrix(int m, int n, double s)
        {
            _m = m;
            _n = n;
            _a = ReturnRectangularDoubleArray(m, n);
            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    _a[i][j] = s;
                }
            }
        }

        // ctor - Construct a matrix from a 2-D array.
        public Matrix(double[][] a)
        {
            _m = a.Length;
            _n = a[0].Length;
            for (int i = 0; i < _m; i++)
            {
                if (a[i].Length != _n)
                {
                    throw new ArgumentException("All rows must have the same length.");
                }
            }
            _a = a;
        }


        // ctor - Construct a matrix quickly without checking arguments
        public Matrix(double[][] a, int m, int n)
        {
            _a = a;
            _m = m;
            _n = n;
        }


        // ctor - Construct a matrix from a one-dimensional packed array of doubles, packed by columns (ala            Fortran).
        public Matrix(double[] vals, int mm)
        {
            _m = mm;
            _n = (_m != 0 ? vals.Length / _m : 0);
            if (_m * _n != vals.Length)
            {
                throw new ArgumentException(
                    "Array length must be a multiple of m.");
            }
            _a = ReturnRectangularDoubleArray(_m, _n);
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    _a[i][j] = vals[i + j * _m];
                }
            }
        } 
        #endregion

        #region ICloneable Members

        public object Clone()
        {
            return Copy();
        }

        #endregion

        #region ISerializable Members

        public void GetObjectData(SerializationInfo info, StreamingContext context)
        {
            info.AddValue("MatrixArray", MatrixArray);
            info.AddValue("RowDimension", RowDimension);
            info.AddValue("ColumnDimension", ColumnDimension);


        }

        #endregion

        private void CheckMatrixDimensions(Matrix matrix)
        {
            if (matrix == null) throw new ArgumentNullException("matrix");
            if (matrix._m != _m || matrix._n != _n)
            {
                throw new ArgumentException("Matrix dimensions must agree.");
            }
        }

        public static Matrix CreateDiagonalMatrix(List<Double> values)
        {
            var m = new Matrix(values.Count, values.Count, 0);
            for (int i = 0; i < values.Count; i++)
            {
                m[i, i] = values[i];
            }
            return m;
        }

        internal static double[][] ReturnRectangularDoubleArray(int size1, int size2)
        {
            var retval = new double[size1][];
            for (int array1 = 0; array1 < size1; array1++)
            {
                retval[array1] = new double[size2];
            }
            return retval;
        }

        // Construct a matrix from a copy of a 2-D array.
        public static Matrix ConstructWithCopy(double[][] a)
        {
            int m = a.Length;
            int n = a[0].Length;
            var x = new Matrix(m, n);
            double[][] c = x.MatrixArray;
            for (int i = 0; i < m; i++)
            {
                if (a[i].Length != n)
                {
                    throw new ArgumentException( "All rows must have the same length.");
                }
                for (int j = 0; j < n; j++)
                {
                    c[i][j] = a[i][j];
                }
            }
            return x;
        }

        // Make a deep copy of a matrix
        public Matrix Copy()
        {
            var x = new Matrix(_m, _n);
            double[][] c = x.MatrixArray;
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    c[i][j] = _a[i][j];
                }
            }
            return x;
        }

        //Clone the Matrix object - return 2D copy of array of matrix elements.
        public double[][] GetArrayCopy()
        {
            double[][] c = ReturnRectangularDoubleArray(_m, _n);
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    c[i][j] = _a[i][j];
                }
            }
            return c;
        }

        // Make a one-dimensional column packed copy of the internal array. return Matrix elements packed in a one-dimensional array by columns.
	    public double[] GetColumnPackedCopy()
        {
            var vals = new double[_m*_n];
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    vals[i + j*_m] = _a[i][j];
                }
            }
            return vals;
        }

       
	    //Make a one-dimensional row packed copy of the internal array - return Matrix elements packed in a one-dimensional array by rows.
	    public double[] GetRowPackedCopy()
        {
            var vals = new double[_m*_n];
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    vals[i*_n + j] = _a[i][j];
                }
            }
            return vals;
        }

        // Get / Set a single element.
        public double this[int i, int j]
        {
            get {return _a[i][j];}
            set { _a[i][j] = value; }
        }

        // TODO: remove GetMatrix()
        public double[][] GetMatrix()
        {
            return _a; 
        }

        // Get a submatrix, given initial / final row & column indices
	    public Matrix GetMatrix(int i0, int i1, int j0, int j1)
        {
            var x = new Matrix(i1 - i0 + 1, j1 - j0 + 1);
            double[][] b = x.MatrixArray;
            for (int i = i0; i <= i1; i++)
            {
                for (int j = j0; j <= j1; j++)
                {
                    b[i - i0][j - j0] = _a[i][j];
                }
            }
            return x;
        }

        // Get a submatrix.
        public Matrix GetMatrix(int[] r, int[] c)
        {
            var x = new Matrix(r.Length, c.Length);
            double[][] b = x.MatrixArray;
            
            for (int i = 0; i < r.Length; i++)
            {
                for (int j = 0; j < c.Length; j++)
                {
                    b[i][j] = _a[r[i]][c[j]];
                }
            }
           
            return x;
        }

        // Get a submatrix: Initial row index, Final row index, Array of column indices.
        public Matrix GetMatrix(int i0, int i1, int[] c)
        {
            var x = new Matrix(i1 - i0 + 1, c.Length);
            double[][] b = x.MatrixArray;
            
            for (int i = i0; i <= i1; i++)
            {
                for (int j = 0; j < c.Length; j++)
                {
                    b[i - i0][j] = _a[i][c[j]];
                }
            }
            
            return x;
        }

        
	    // Get a submatrix  - Array of row indices, Initial column index,  Final column index
        public Matrix GetMatrix(int[] r, int j0, int j1)
        {
            var x = new Matrix(r.Length, j1 - j0 + 1);
            double[][] b = x.MatrixArray;
            try
            {
                for (int i = 0; i < r.Length; i++)
                {
                    for (int j = j0; j <= j1; j++)
                    {
                        b[i][j - j0] = _a[r[i]][j];
                    }
                }
            }
            catch (IndexOutOfRangeException)
            {
                throw new IndexOutOfRangeException("Submatrix indices");
            }
            return x;
        }

        // Set a submatrix - Initial row index, Final row index, Initial column index, Final column index
	    public void SetMatrix(int i0, int i1, int j0, int j1, Matrix x)
        {
            for (int i = i0; i <= i1; i++)
            {
                for (int j = j0; j <= j1; j++)
                {
                    _a[i][j] = x[i - i0, j - j0];
                }
            }
           
        }

       // Set a submatrix -  Array of row indices.,  Array of column indices.
	    public void SetMatrix(int[] r, int[] c, Matrix x)
        {
            for (int i = 0; i < r.Length; i++)
            {
                for (int j = 0; j < c.Length; j++)
                {
                    _a[r[i]][c[j]] = x[i, j];
                }
            }
        }

     
        // Set a submatrix -  Array of row indices,  Initial column index, Final column index
        public void SetMatrix(int[] r, int j0, int j1, Matrix x)
        {
            for (int i = 0; i < r.Length; i++)
            {
                for (int j = j0; j <= j1; j++)
                {
                    _a[r[i]][j] = x[i, j - j0];
                }
            }
        }

       
        // Set a submatrix -  Initial row index, Final row index, Array of column indices,  
        public void SetMatrix(int i0, int i1, int[] c, Matrix x)
        {
            for (int i = i0; i <= i1; i++)
            {
                for (int j = 0; j < c.Length; j++)
                {
                    _a[i][c[j]] = x[i - i0, j];
                }
            }
        }

        //  Matrix transpose.
        public Matrix Transpose()
        {
            var x = new Matrix(_n, _m);
            double[][] c = x.MatrixArray;
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    c[j][i] = _a[i][j];
                }
            }
            return x;
        }

      
	    // One norm - return maximum column sum.
        public double Norm1()
        {
            double f = 0;
            for (int j = 0; j < _n; j++)
            {
                double s = 0;
                for (int i = 0; i < _m; i++)
                {
                    s += Math.Abs(_a[i][j]);
                }
                f = Math.Max(f, s);
            }
            return f;
        }

        /**
	 * Two norm - @return maximum singular value.
	 */

        // public double norm2 () {
        // return (new SingularValueDecomposition(this).norm2());
        // }

        // Infinity norm - return maximum row sum.
        public double NormInf()
        {
            double f = 0;
            for (int i = 0; i < _m; i++)
            {
                double s = 0;
                for (int j = 0; j < _n; j++)
                {
                    s += Math.Abs(_a[i][j]);
                }
                f = Math.Max(f, s);
            }
            return f;
        }

        /**
	 * Frobenius norm
	 * 
	 * @return sqrt of sum of squares of all elements.
	 */

        // public double normF () {
        // double f = 0;
        // for (int i = 0; i < m; i++) {
        // for (int j = 0; j < n; j++) {
        // f = Maths.hypot(f,A[i][j]);
        // }
        // }
        // return f;
        // }
     

        // Unary minus
        public Matrix Uminus()
        {
            var x = new Matrix(_m, _n);
            double[][] c = x.MatrixArray;
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    c[i][j] = -_a[i][j];
                }
            }
            return x;
        }

        /**
	 * C = A + B
	 * 
	 * @param B
	 *            another matrix
	 * @return A + B
	 */
        
        public Matrix Plus(Matrix b)
        {
            CheckMatrixDimensions(b);
            var x = new Matrix(_m, _n);
            double[][] c = x.MatrixArray;
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    c[i][j] = _a[i][j] + b._a[i][j];
                }
            }
            return x;
        }

      
	
	    // @return A + B
        public Matrix PlusEquals(Matrix b)
        {
            CheckMatrixDimensions(b);
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    _a[i][j] = _a[i][j] + b._a[i][j];
                }
            }
            return this;
        }

        // @return A - B
        public Matrix Minus(Matrix b)
        {
            CheckMatrixDimensions(b);
            var x = new Matrix(_m, _n);
            double[][] c = x.MatrixArray;
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    c[i][j] = _a[i][j] - b._a[i][j];
                }
            }
            return x;
        }

        // return A - B
        public Matrix MinusEquals(Matrix b)
        {
            CheckMatrixDimensions(b);
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    _a[i][j] = _a[i][j] - b._a[i][j];
                }
            }
            return this;
        }

        // return A * B
        public Matrix ArrayTimes(Matrix b)
        {
            CheckMatrixDimensions(b);
            var x = new Matrix(_m, _n);
            double[][] c = x.MatrixArray;
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    c[i][j] = _a[i][j]*b._a[i][j];
                }
            }
            return x;
        }

        // return A * B
        public Matrix ArrayTimesEquals(Matrix b)
        {
            CheckMatrixDimensions(b);
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    _a[i][j] = _a[i][j]*b._a[i][j];
                }
            }
            return this;
        }

        // @return A / B
        public Matrix ArrayRightDivide(Matrix b)
        {
            CheckMatrixDimensions(b);
            var x = new Matrix(_m, _n);
            double[][] c = x.MatrixArray;
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    c[i][j] = _a[i][j]/b._a[i][j];
                }
            }
            return x;
        }

        // return A / B
        public Matrix ArrayRightDivideEquals(Matrix b)
        {
            CheckMatrixDimensions(b);
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    _a[i][j] = _a[i][j]/b._a[i][j];
                }
            }
            return this;
        }

        // return A / B
        public Matrix ArrayLeftDivide(Matrix b)
        {
            CheckMatrixDimensions(b);
            var x = new Matrix(_m, _n);
            double[][] c = x.MatrixArray;
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    c[i][j] = b._a[i][j]/_a[i][j];
                }
            }
            return x;
        }

        // @return A / B
        public Matrix ArrayLeftDivideEquals(Matrix b)
        {
            CheckMatrixDimensions(b);
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    _a[i][j] = b._a[i][j]/_a[i][j];
                }
            }
            return this;
        }

        // return s * A 
        public Matrix Times(double s)
        {
            var x = new Matrix(_m, _n);
            double[][] c = x.MatrixArray;
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    c[i][j] = s*_a[i][j];
                }
            }
            return x;
        }

        // @return s * A 
        public Matrix TimesEquals(double s)
        {
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    _a[i][j] = s*_a[i][j];
                }
            }
            return this;
        }

        // @return A * B
        public Matrix Times(Matrix b)
        {
            if (b._m != _n)
            {
                throw new ArgumentException("Matrix inner dimensions must agree.");
            }
            var x = new Matrix(_m, b._n);
            double[][] c = x.MatrixArray;
            var bcolj = new double[_n];
            for (int j = 0; j < b._n; j++)
            {
                for (int k = 0; k < _n; k++)
                {
                    bcolj[k] = b._a[k][j];
                }
                for (int i = 0; i < _m; i++)
                {
                    double[] arowi = _a[i];
                    double s = 0;
                    for (int k = 0; k < _n; k++)
                    {
                        s += arowi[k]*bcolj[k];
                    }
                    c[i][j] = s;
                }
            }
            return x;
        }

        /**
	 * LU Decomposition
	 * 
	 * @return LUDecomposition
	 * @see LUDecomposition
	 */

        public LuDecomposition Lu()
        {
            return new LuDecomposition(this);
        }

        // /** QR Decomposition
        // @return QRDecomposition
        // @see QRDecomposition
        // */
        //
        // public QRDecomposition qr () {
        // return new QRDecomposition(this);
        // }
        //
        // /** Cholesky Decomposition
        // @return CholeskyDecomposition
        // @see CholeskyDecomposition
        // */
        //
        // public CholeskyDecomposition chol () {
        // return new CholeskyDecomposition(this);
        // }
        //
        // /** Singular Value Decomposition
        // @return SingularValueDecomposition
        // @see SingularValueDecomposition
        // */
        //
        // public SingularValueDecomposition svd () {
        // return new SingularValueDecomposition(this);
        // }
        //
        // /** Eigenvalue Decomposition
        // @return EigenvalueDecomposition
        // @see EigenvalueDecomposition
        // */
        //
        // public EigenvalueDecomposition eig () {
        // return new EigenvalueDecomposition(this);
        // }

        /**
	 * Solve A*X = B
	 * 
	 * @param B
	 *            right hand side
	 * @return solution if A is square, least squares solution otherwise
	 */

        public Matrix Solve(Matrix b)
        {
            // assumed m == n
            return new LuDecomposition(this).Solve(b);
        }

       // Solve X*A = B - return solution  if A is square, least squares solution otherwise.
	    public Matrix SolveTranspose(Matrix b)
        {
            return Transpose().Solve(b.Transpose());
        }

        // Matrix inverse 
        public Matrix Inverse()
        {
            return Solve(Identity(_m, _m));
        }

     
	    // Matrix determinant
        public double Det()
        {
            return new LuDecomposition(this).Det();
        }

        // Matrix rank - return effective numerical rank, obtained from SVD.
	    // public int rank () {
        // return new SingularValueDecomposition(this).rank();
        // }
        
        // Matrix condition (2 norm) - return ratio of largest to smallest singular value.
        // public double cond () {
        // return new SingularValueDecomposition(this).cond();
        // }

        // Matrix trace - return sum of the diagonal elements.
        public double Trace()
        {
            double t = 0;
            for (int i = 0; i < Math.Min(_m, _n); i++)
            {
                t += _a[i][i];
            }
            return t;
        }

       // Generate matrix with random elements - return An m-by-n matrix with uniformly distributed random elements.
	    public static Matrix Random(int m, int n)
        {
            var a = new Matrix(m, n);
            double[][] x = a.MatrixArray;
            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    x[i][j] = new Random(1).NextDouble();
                }
            }
            return a;
        }

     
	    // Generate identity matrix
        public static Matrix Identity(int m, int n)
        {
            var a = new Matrix(m, n);
            double[][] x = a.MatrixArray;
            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    x[i][j] = (i == j ? 1.0 : 0.0);
                }
            }
            return a;
        }

        /**
	 * Print the matrix to stdout. Line the elements up in columns with a
	 * Fortran-like 'Fw.d' style format.
	 * 
	 * @param w
	 *            Column width.
	 * @param d
	 *            Number of digits after the decimal.
	 */

        // TODO: Matrix.print - PrintWriter

        //public void print(int w, int d) {
        //    print(new PrintWriter(System.out, true), w, d);
        //}

        /**
	 * Print the matrix to the output stream. Line the elements up in columns
	 * with a Fortran-like 'Fw.d' style format.
	 * 
	 * @param output
	 *            Output stream.
	 * @param w
	 *            Column width.
	 * @param d
	 *            Number of digits after the decimal.
	 */
        // TODO: Matrix.print - PrintWriter, DecimalFormat
        //public void print(PrintWriter output, int w, int d) {
        //    DecimalFormat format = new DecimalFormat();
        //    format.setDecimalFormatSymbols(new DecimalFormatSymbols(Locale.US));
        //    format.setMinimumIntegerDigits(1);
        //    format.setMaximumFractionDigits(d);
        //    format.setMinimumFractionDigits(d);
        //    format.setGroupingUsed(false);
        //    print(output, format, w + 2);
        //}

        /**
	 * Print the matrix to stdout. Line the elements up in columns. Use the
	 * format object, and right justify within columns of width characters. Note
	 * that is the matrix is to be read back in, you probably will want to use a
	 * NumberFormat that is set to US Locale.
	 * 
	 * @param format
	 *            A Formatting object for individual elements.
	 * @param width
	 *            Field width for each column.
	 * @see java.text.DecimalFormat#setDecimalFormatSymbols
	 */

        // TODO: Matrix.print - PrintWriter, NumberFormat
        //public void print(NumberFormat format, int width) {
        //    print(new PrintWriter(System.out, true), format, width);
        //}

        // DecimalFormat is a little disappointing coming from Fortran or C's
        // printf.
        // Since it doesn't pad on the left, the elements will come out different
        // widths. Consequently, we'll pass the desired column width in as an
        // argument and do the extra padding ourselves.

        /**
	 * Print the matrix to the output stream. Line the elements up in columns.
	 * Use the format object, and right justify within columns of width
	 * characters. Note that is the matrix is to be read back in, you probably
	 * will want to use a NumberFormat that is set to US Locale.
	 * 
	 * @param output
	 *            the output stream.
	 * @param format
	 *            A formatting object to format the matrix elements
	 * @param width
	 *            Column width.
	 * @see java.text.DecimalFormat#setDecimalFormatSymbols
	 */
        // TODO: Matrix.print - PrintWriter, NumberFormat
        //public void print(PrintWriter output, NumberFormat format, int width) {
        //    output.println(); // start on new line.
        //    for (int i = 0; i < m; i++) {
        //        for (int j = 0; j < n; j++) {
        //            String s = format.format(A[i][j]); // format the number
        //            int padding = Math.Max(1, width - s.Length()); // At _least_ 1
        //            // space
        //            for (int k = 0; k < padding; k++)
        //                output.print(' ');
        //            output.print(s);
        //        }
        //        output.println();
        //    }
        //    output.println(); // end with blank line.
        //}

        /**
	 * Read a matrix from a stream. The format is the same the print method, so
	 * printed matrices can be read back in (provided they were printed using US
	 * Locale). Elements are separated by whitespace, all the elements for each
	 * row appear on a single line, the last row is followed by a blank line.
	 * 
	 * @param input
	 *            the input stream.
	 */
        // TODO: Matrix.read --> deserializer
        //@SuppressWarnings("unchecked")
        //public static Matrix read(BufferedReader input) throws java.io.IOException {
        //    StreamTokenizer tokenizer = new StreamTokenizer(input);

        //    // Although StreamTokenizer will parse numbers, it doesn't recognize
        //    // scientific notation (E or D); however, Double.valueOf does.
        //    // The strategy here is to disable StreamTokenizer's number parsing.
        //    // We'll only get whitespace delimited words, EOL's and EOF's.
        //    // These words should all be numbers, for Double.valueOf to parse.

        //    tokenizer.resetSyntax();
        //    tokenizer.wordChars(0, 255);
        //    tokenizer.whitespaceChars(0, ' ');
        //    tokenizer.eolIsSignificant(true);
        //    java.util.Vector v = new java.util.Vector();

        //    // Ignore initial empty lines
        //    while (tokenizer.nextToken() == StreamTokenizer.TT_EOL)
        //        ;
        //    if (tokenizer.ttype == StreamTokenizer.TT_EOF)
        //        throw new java.io.IOException("Unexpected EOF on matrix read.");
        //    do {
        //        v.addElement(Double.valueOf(tokenizer.sval)); // Read & store 1st
        //        // row.
        //    } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);

        //    int n = v.Count; // Now we've got the number of columns!
        //    double row[] = new double[n];
        //    for (int j = 0; j < n; j++)
        //        // extract the elements of the 1st row.
        //        row[j] = ((Double) v.elementAt(j)).doubleValue();
        //    v.removeAllElements();
        //    v.addElement(row); // Start storing rows instead of columns.
        //    while (tokenizer.nextToken() == StreamTokenizer.TT_WORD) {
        //        // While non-empty lines
        //        v.addElement(row = new double[n]);
        //        int j = 0;
        //        do {
        //            if (j >= n)
        //                throw new java.io.IOException("Row " + v.Count
        //                        + " is too long.");
        //            row[j++] = Double.valueOf(tokenizer.sval).doubleValue();
        //        } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
        //        if (j < n)
        //            throw new java.io.IOException("Row " + v.Count
        //                    + " is too short.");
        //    }
        //    int m = v.Count; // Now we've got the number of rows.
        //    double[][] A = new double[m][];
        //    v.copyInto(A); // copy the rows out of the vector
        //    return new Matrix(A);
        //}

        public override String ToString()
        {
            var buf = new StringBuilder();
            for (int i = 0; i < RowDimension; i++)
            {
                for (int j = 0; j < ColumnDimension; j++)
                {
                    buf.Append(this[i, j]);
                    buf.Append(" ");
                }
                buf.Append("\n");
            }

            return buf.ToString();
        }
    }
}